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ABSTRACT 



Aims. Over the past several years, numerous examples of X-ray cavities coincident with radio sources have been observed in so- 
called "cool core" clusters of galaxies. Motivated by these observations, we explore the evolution and the effect of cavities on a 
cooling intracluster medium (ICM) numerically, adding relevant physics step by step. 

Methods. In this paper we present a first set of hydrodynamical, high resolution (1024^^ effective grid elements), three-dimensional 
simulations, together with two-dimensional test cases. The simulations follow the evolution of radio cavities, modeled as bubbles 
filled by relativistic plasma, in the cluster atmosphere while the ICM is subject to cooling. 

Results. We find that the bubble rise retards the development of a cooling flow by inducing motions in the ICM which repeatedly 
displace the material in the core. Even bubbles initially set significantly far from the cluster center affect the cooling flow, although 
much later than the beginning of the simulation. The effect is, however, modest: the cooling time is increased by at most only 25%. 
As expected, the overall evolution of pure hydrodynamic bubbles is at odds with observations, showing that some additional physics 
has to be considered in order to match the data. 

Key words. Galaxies; clusters; general - cooling flows - methods; numerical 



1. Introduction 

Approaching the centers of rich clusters of galaxies, one often 
finds a remarkable increase in the X-ray luminosity together 
with a marked decrease in the temperature. Under the condi- 
tions present in these clusters, the radiative cooling time of 
the intracluster medium (ICM) becomes shorter than the clus- 
ter lifetime and a "cooling flow" is expected to occur, i.e. a 
subson i c flow of cooling gas toward the cluster center (see 
iFabianl 11994 for a review). The absence of observations of 
a reliable repository for the cold gas left some doubts regard- 
ing this picture, which were confirmed when a major revi- 
sion was forced by spectral measuremen ts of the cores of the 
cluste r s. High-resolu t ion XMM/Newton (Peterson et all |2001L 



l2003t iTamura et al.l |200U 'Kaastra et al 



lower resolu t ion A SCA ^akishima et al.i 12001 



2004 



spectra and 
) and Chandra 



jPavid et al.L 12001') spectra show that only a small fraction of 
the gas is cooling below 1-2 keV, whereas its temperature was 
expected to be as low as 10^ K (0.1 keV). Hence it is now said 
that these clusters host a "cool core" instead of a cooling flow, 
and therefore in the rest of the paper we will refer to them as 
cool core cluster 

Maintaining gas at keV temperatures for a period several 
times longer than its cooling time requires one or more heat- 
ing mechanisms. Candidate heating mechanisms include ther- 
mal c ond uction by electron s in the ICM (Narayan & Medvede^ 
1200 It IZakamska & Naravan , 2003; Voigt & Fabian , 2004), re- 
connection o f magnetic fields (.Soker & Sarazin. .1990) , turbu- 
lent mixing dKim & Naravanl l2003h . and the injection of hot 
plasma by the active galactic nucleus (AGN) hosted by the cen- 
tral galaxy. In this paper we will focus on this last possibility, but 



we do not exclude the possibility that more than one mechanism 
is at work. 

Recent observations suggest that almost all c ool co re clus- 
ters host radio sources in their centers (EileJ, l2004 . Often 
these sources are rather faint, and it is not clear if a correla- 
tion exists between their radio power and the strength of the ex- 
pected cooUng flo\y dVoigt & Fabian l l2004lKaastra et al.Ll2004 
iBirzan et al.L |2004) . However, about 71% of the cD galaxies in 
cool core clusters are r adio -loud com pared to only 23 % of non- 
cool core cluster cDs (Burns, 1990; Ball etaUll993l) . and this 
suggests a connection between AGN activity and the presence 
of a cool core, at least in clusters that host a cD. 

The shapes of radio sources are determined by the interaction 
of the jets with their surroundings. The central radio sources in 
cool core clusters mostly show a disturbed morphology: even 
when they host a radio-loud core, coUimated jets exist only on 
kpc scales or below. After this the plasma flow continues in a less 
collimated manner into the ICM, and the lobes take the shape 
of lar ge plumes or wide tails (see e.g. lEilekl |2004 iJetha et akl 
l2005h . A remarkable exception is Cygnus A, a typical Fanaroff- 
Riley II (FRII) source, whose jets travel straight for about 100 
kp c and terminate in hot spo ts at the far ends of the lobes (see 
e.g lBardiel & Arnaudl [19961) . 

The interaction of a single jet or a pair of jets with 
the ICM has been studied extensively both theoretically 
(Scheuer, 1974; Blandford & Rees, 1974; Smith etal., 1983|; 
Begelman & Cioffi. .1989: .Heinz et al.. 1998; Alexander, 200% 
and through nume r ical simulations (see e .^. Clarke et al., 1993 

2002t iKrausei i2003t 



'Rizza et al.', '2000; 'Reynolds et alj |2001 



JCrause & Camenzind, 2003; Zan ni et al 



2003L and references 
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therein). These simulations, however, mostly apply only to FRII 
sources and to the phase of AGN activity. This work sug- 



2 



Gardini: Buoyant Bubbles in a Cooling ICM 



gests that jets inflate a cocoon with radio plasma which en- 
compasses both the jets and the AGN, and as the radio source 
evolves, this cocoon elongates in the jet's direction and as- 
sumes a cigar-like shape. The long-term evolution of the ra- 
dio sources has received less attention; on the one hand, it 
has been argued that they should continue to expand sub- 
sonically while the cocoon remains overpressured with re- 
spect to the surrounding medium. On the other hand, it was 
pointed out that the cocoon evolution would instead be driven 
by buoyancy when its expan sion speed became co mparable 
to the rising buoyancy speed jChurazov et al.L l2000h . A more 
detailed picture has been reveal e d by some recent numeri- 
cal studies ( Revn olds eTaLl l2002t |B asson & A lexandeR l2003t 
lOmma et all 12004; O mma & Binnevi 12004V Zanni et al. , 20051): 
after the jets are switched ofi^, they show that lateral instabilities 
bisect the cocoon in its central region, while its remaining parts 
start to rise in the cluster atmosphere, due to residual momentum 
and to buoyancy. 

Depressions of the X-ray flux coincident with AGN radio 
lobes have been observed in cool core clusters and also in some 
groups and galaxies by ROSAT and Chandra telescopes. They 
are usually interpreted as cavities in the ICM, produced and filled 
by the radio-emitting plasma injected by the central AGN as de- 
scr ibed above. A systematic study of these systems is provided 
bv lBirzan et al.l (|2004). Usually, cavities occur in pairs, and the 
cavities in the same pair are in opposite directions and at about 
the same distance with respect to the cluster center. In general, 
they also have similar dimensions, measuring up to some tens 
of kpc across, and are located some tens of kpc away from the 
cluster center However, a pair of huge cavities, roughly 200 
kpc in diameter, has been detected by McNamara et al. (2005.) 
in MS0735.6H-7421. 

Cavities are sometimes surrounded by bright rims of denser 
and cooler gas (see e.g. A2052: Blanton et al., 2001), which pre- 
sumably has been displaced, uplifted or entrained by the hot 
plasma during the formation of the cavity or its subsequent mo- 
tion. In some cases weak shocks are observed surrounding the 
cavities in clusters (Perseus: Fabian et al., 2006, 2003a; Hydra 
A: Nulsen et al., 2005b; MS0735.6+7421: McNamara et al., 
2005) and in M87 (iForman et al.Ll2005l) . b ut so far strong sho cks 
have been detected only in Centaurus A jKraft et all l2003l) . A 
weak shock has also been report ed around the radio source in 
Hercules A dNulsen et al.L l2005 a). which apparently does not 
host c avities, and spec t ral ev idence of a shock is reported in 
A478 dSanderson et all |2005|) . A bow shoc k in front of radio 
lobes may have been detected i n Cygnus A (ISmifli et al.Ll2002t 
iBalucinska-Church et al.L l2005l) . In the case of a weak shock, it 
can be assumed that the radio-emitting plasma in the cavities is 
almost in pressure equilibrium with the surrounding ICM. Thus, 
due to their low density, the radio lobes will rise by buoyancy 
in the cluster atmosphere (iGuU & Northoverlll973h . similarly to 
the remnants of the cocoon in the models described above. 

Chandra X-ray images also show cavities that are not co- 
incident with bright rad io lobes, fo r example in Abell 2597 
(iMcNamara etaLl 120011). Per seus (iFabian et al.L l2000l) . and 
Abell 4059 dHeinz et al. , 1200 2*). These structures are referred to 
as "ghost cavities". The relativistic electrons in the radio lobes 
are expected to lose enough energy via synchrotron emission to 
become invisible in the radio after 50 to 100 Myr. If this interpre- 
tation is correct, the ghost cavities are buoyantly risi ng relics of a 
radio outburst that ended at least 50 to 100 Myr ago (iSoker et al.L 
I2002h . 

The commonly observed "cold fronts" discovered by 
Chandra can also provide hints regarding the effect of AGN 



on the intracluster medium. Cold fronts are sharp disconti- 
nuities in X-ray surface brightness marking boundaries be- 
tween hotter and colder masses of gas in the ICM. The den- 
sity and temperature jumps inferred for cold fronts through 
deprojection analysis are consistent with continuous pressure 
changes acro ss the fronts, showing th at the discontinuities are 
not shocks ( Markevitch et al.L l2000h . While cold fronts are 
usually explained as remnants of past merger events, they 
have been observed als o in apparently regular clusters such 
as RXJ 1720.1+ 2638 d Mazzottaet a l.l. l200 l ab, Abell 179 5 
dMarkevitch e t al.L l200ll). 2A 0335-H09 6 (Mazz otta et aill2003l) . 
and MS 1455.0+2232 ('Maz zotta eTan.l2001 b). The gas sloshing 
induced by th e motion of a pre-exist ing cavity may explain these 
observations dMazzotta et al.Ll2003h . 

The Ha filaments observed in the core of the Perseus cluster 
also appear to be related to the upward motion of the cavities. 
In particular, some filaments seem drawn up in a laminar flow 
in the wake of a ghost cavity, also tracing well-defined arcs that 
resemble a circulation flow. This suggests that t he medium is 
not turbulent and has a non-negligible viscosity dFabian et all 
2003b). 

Studies of the heating of the ICM through energy 
injection by AGN have also been performed using an- 
alytical and semi-anal yti cal m o dels (Beaelman, 2001; 
Ruszkowski & Begelman, 2002; Bohringer et al., 2002; 
Chur azoveLa l . , 2002; Kaiser & B inney, 2003; Briiggen, 2Q03al ; 
Ma thews et al. . 2003, 2004, 20061: iRoychowdhury et all l2004 . 
which, however, describe these systems in very general terms. 
This is due to the fact that the dynamics of the hot plasma 
in the ICM is a very complex problem, owing to its intrinsic 
three-dimensionality, and to the chaotic interactions between 
plasma and ICM. W hile analytical models of buoyant bubbles 
have bee n produc ed ("Gull & Northo veJ. 119731: IChurazov et all 
2000; So ker et all 2002; Kaiser et al.r i2005l) . the main tool of 
investigation consists of numerical simulations. The interaction 
of the radio sources with the ICM has been simulated so far in 
different ways. Their formation is often m odeled by injecting 
energy or hot plasma directly into the ICM ('Ouilis et al., 20011; 
Briahenti & Mathews, 2002a b, 2003; Bru2aen & Kaiser, 2Q(A 
Briiggen et al., 2002, 2005; Briiggen, 2003b; Ruszkowsk i et all 
200 4a.k .DaUa Vecchia et al.. .2004: Jones & De Young. |2005[ 
Ver naleo & Reynolds! l2006h . rather than reproducing detailed 
AGN jets as discussed above. This approach enables one to 
study repeated bursts of activity, whereas jet simulations at high 
resolution are usually unable to cover long evolutionary times. 
In a different approach, the phase of buoyancy is reproduced 
by modelUng the radio sources as bubbles o f hot plasma 
(Churazov et al., 2001; Bruggen & Kaisei', '2001'; S axton et all 
2001; Rob inso n et al., 2004; Reynolds et al., 2005), tiius ne- 
glecting the phase of formation of cavities, and starting from a 
configuration in pressure equilibrium similar to observations. 
A review of simu lations of hot plasma in clusters is given in 
iGardini & Ricked d2004 . 

We aim to investigate the properties of the ICM, the dynam- 
ical evolution of the radio cavities and their impact on the envi- 
ronment, by means of simulations of increasing complexity and 
realism. In this paper we show the results of a first set of three- 
dimensional (3D) hydrodynamical simulations in which cavities 
are modeled as buoyant bubbles in the atmosphere of a cool core 
cluster, together with two-dimensional (2D) test cases. With re- 
spect to previous similar numerical experiments, we achieve a 
rather high effective numerical resolution (1024^ grid elements), 
and account at the same time for the cooling of the gas. Similar 
simulations have been performed in three dimensions only by 
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[Reynolds et alj (l2005h . who focus on the effect of ICM viscosity 
on the bubble evolution and therefore neglect cooling. This work 
deals both with the dynamics of the bubbles and their effect on 
the cooling of the ICM. We also stress the reliability and robust- 
ness of the results, in order to use them as reference for future 
simulations. 

The paper is organized as follows: in §2 we describe the char- 
acteristics of the simulations and of the code. In §3 we describe 
the evolution of the simulations, while the effects of bubbles on 
the cool core are discussed in §4. Finally, in §5 we summarize 
our results and discuss future developments of this work. 




2. Characteristics of simulations 



100 
(kpc) 



We defined the cluster parameters in order to resemble the ob- 
served characteristics of Abell 2390, a massive galaxy cluster 
with a significant c ool core. We model the gravita tional potential 
according to the Navarro. Frenk & Wliitd ( 1997b density distri- 
bution for dark-matter halos. 



p(r) = P.S 



(1) 



and keep it fixed during each run. We set the scale radius rs to 
520 kpc, while the scale density ps = 7.07 x 10"^^ g/cm^ is de- 
rived by imposing the requirement that the virial mass (M2(K)) of 
the cluster be equal to 1 .7 x 10''' Mq, after assuming the Hubble 
parameter h - 0.7 and redshift z = 0; the corresponding virial 
radius is then r2oo ~ 1.72/i"' Mpc. 

The initial gas temperature profile is shown in Figure 1 ; it is 
constructed a s an analytic adaptatio n of the profile of Abell 2390 
proposed by Zakamska & NarayanI (120031 which in turn repro- 
duces by eye the temperature profile obtained by deprojection 
from Chandra observations by lAllen et al.l (1200 lb . The expres- 
sion for T(r) is 

T{r) = - (To,, - r,„)exp(-r2/(2(r(r)2)))x 
exp Ulog(r) - log(r,)) - ^//2(log(r) - log(r,))2 + O.OO4]] (2) 



Fig. 1. The temperature profile of the simulated cluster The 
temperature rises from about 4 keV at the center up to 11 
keV at the scale radius, before decre asing in the o utskir ts. The 
profile is derived from the data of lAllen et aP (l200lh as in 
Zakams ka & NarayanI (|2003|) for Abell 2390. Beyond the scale 
ra dius the profile goes as T <x r ~° '*, as suggested by the results 
of jPe Grandi & Molendil (l2002h . 



which requires that the temperature rise from r,„ = 4 keV at the 
center up to Tout = 11 keV at the scale radius. Beyond the scale 
radius we assume that the gas temperature drops like T oc r^,// - 
-0.4, as suggested by the results of iDe Grandi & Molendil 
(I2OO2I) . The definition of cr is 



cr(r) — ar + b 
with 



(3 V-21og0.5j 



/(r., - rave) 



and 



b^rA--a 



(3) 



(4) 



(5) 



where the radius rave = 80 kpc encompasses the part of the clus- 
ter where T < (Tout + Tu,)/2. 

The gas density profile is obtained by solving the equation 
of hydrostatic equilibrium with the condition that the gas frac- 
tion inside the virial radius be ft = 0.15, approaching the cosmic 



value (see e.g. lSpergel et al.Ll2003l) . In Figure 2 we show the pro- 
file of the cooling time fcooi = -(d In T/dt y^ for the initial con- 
ditions according to the cooling function of lSutherland & Dopital 
(Il993h used in the code. It is worth noting that the region of the 
cluster where fcooi ^ 1 Gyr spans only r < 12 kpc. The adia- 
batic index of the ICM is set to y = 5/3, assuming a monatomic 
perfect gas. We tested the equilibrium of the initial configuration 
of the ICM by running a 2D simulation without cooling or bub- 
bles. After f = 500 Myr, readjustments in the physical quantities 
(density, pressure, etc.) amounted at most to 1%. 

The rising radio lobes are simulated as a pair of spherical 
bubbles in the ICM, placed symmetrically with respect to the 
cluster center. Inside the bubbles, we set the gas density and pres- 
sure respectively to 1/100 of and the same value as the ICM's 
quantities at the same radius; we also impose y - 4/3 in the 
bubble material because the plasma is relativistic. We did not at- 
tempt to include the ICM displaced from these cavities during 
the phase of bubble formation. In the following we will refer to 
the two gas components which had y = 5/3 and 7 - 4/3 respec- 
tively as ICM and bubble plasma. 

The different characteristics of the 3D runs are summarized 
in Table 1. The letters CO (cooling only) identify simulations 
performed without bubbles. In simulations with bubbles, the ini- 
tials SB (single bubble pair) are followed by the radius Ri, - 10 
kpc of the bubbles and the distance di, of their centers from the 
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Table 1. Parameters of the 3D simulations. 



Name 




Ax" 








f filial 


CO 


4800 


4688 








500 


C02 


600 


586 








250 


SBlOll 


600 


586 


10 


11 


337 


250 


SB 1020 


600 


586 


10 


20 


281 


250 


SB 1030 


600 


586 


10 


30 


235 


250 



" The box size in kpc 

* The maximum resolution in pc 



The radius of the bubbles in kpc 

The distance of the center of the bubbles from the center of the 
cluster in kpc 

*■ The internal energy of the bubbles in units of 10^^' erg 
^ The time of the final output of the simulation in Myr 



Fig. 2. The value of the cooling time tcooi - -(d In T/dt) ' as a 
function of the radius. 



center of the cluster. The increasing computational cost did not 
permit us to run the simulations for their entire cooling time, but 
they cover a significant fraction of it, up to f = 250 Myr. A set 
of simulations with Rb - 20 kpc halted at earlier stages of evo- 
lution and has not been included in this analysis. All of these 
simulations are performed on a Cartesian grid. We also created a 
number of 2D simulations in cylindrical coordinates to test and 
extend the 3D results. Because of the symmetry of the system, it 
can be noted that the same amount of information extracted by 
a 3D simulation could be obtained through a faster and compu- 
tationally cheaper 2D counterpart. We noticed, however, that the 
two bubbles in these 2D runs evolve in a similar but not fully 
symmetric manner. This is probably due to a greater amplifica- 
tion of rounding errors in this particular system of coordinates 
with respect to the Cartesian grid. The evolution of the bubbles 
and their effect on the medium are nonetheless in agreement with 
the 3D cases until they are available, and this encouraged us to 
trust the 2D results, at least as checks or hints of the overall evo- 
lution of the system. We will not discuss the 2D runs extensively 
but we will report their results when necessary. 

The 3D computational box is centered on the cluster center 
and in general spans 600 kpc; its dimensions are a good compro- 
mise between resolution and the need to avoid spurious bound- 
ary effects. Because we chose reflecting boundary conditions to 
conserve mass and energy, the box size has to be large enough 
that the cooling flow and the bubbles are unaffected. We achieve 
a maximum resolution of ~ 586 pc, comparable to the electron 
mean fr ee path in the c luster core, when magnetic fields are ne- 
glected (lSarazinlll988l) . The CO model, which uses a box size 



of 4800 kpc and a resolution of ~ 4.7 kpc, is a test case for the 
development of the cooling flow in a larger box. 

The cooling function of [Sutherland & Dopital ( 1 19931) is im- 
plemented in the code; we shut off cooUng below a temperature 

of 10"^ K and assume an ICM metallicity of ~ Q.3Za- 

All simulations are performed using FLASH (iFrvxell et al.L 

l2000h . a parallel Eulerian hydrodynamic co de which implements 
the pie cewise-parabolic method (PPM) o f Colella & Woodward I 
(11984 on an adaptive mesh. Grid cells are grouped in blocks, 
which are refined or derefined according to the maximum valu e 
of the second derivative of pressure or density (lLohnerLll987h . 
In 3D simulations the volume is divided into 4^ blocks at the 
coarsest level of refinement, and each block contains 8"^ cells. 
Up to 6 levels of refinements are allowed, with a factor of two 
refinement between levels. 

FLASH manages mixtures of gases using the method of 
IColella & GlazI (Il985h to handle general equations of state. For 
the purposes of this work it is worth noting that each fluid el- 
ement in the simulation is composed of a fraction of both the 
ICM and the bubble plasma, and that all of the physical quan- 
tities, such as density, temperature and pressure, are evaluated 
on the mixture of the two species. In particular to each cell it is 
assigned a weighted average adiabatic index y as 

(r-1) " in - 1) 

where X, is the mass fraction of the /th species and y,- its adiabatic 
index. This left us with two possible problems: on one hand, 
mixing is not properly treated in these simulations because the 
diffusion of species is not explicitly accounted for. On the other 
hand, the temperature in the plasma elements drops quickly as 
they are polluted by the ICM, but the relativistic value y = 4/3 
of their adiabatic index does not change. We address the possible 
effects of these inconsistencies in the following sections. 

3. Evolution of simulations 

3. 1 . The development of the cooling flow 

The simulations named CO and C02 in Table 1 allow the ICM 
to cool, in the absence of bubbles or any source of heating. They 
reproduce the effect of cooling on the ICM starting from the 
equilibrium configuration assumed in the initial conditions. The 
simulations are identical except for the box size and the corre- 
sponding resolution. 

Figure 3 shows the evolution of temperature T, electron den- 
sity n, entropy s - T/n^^^, and pressure P in the central part 
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10 100 10 100 

r(kpc} r(kpc) 



Fig. 3. The evolution of the profiles of some of the ICM quan- 
tities due to cooling, in the absence of cavities. The plot shows 
the evolution of temperature T, electron density n, entropy s - 
T jrp-l^, and pressure P in the central part of the cluster. Profiles 
are obtained from the C02 (solid lines) and CO (dashed lines) 
simulations by averaging over radial bins, assuming a bin inter- 
val of 3 kpc. The time evolution for the C02 simulation is repre- 
sented by the line thickness, which decreases with time. Profiles 
are computed with an interval of 50 Myr from the initial condi- 
tions up to f = 250 Myr The dashed lines represent the profiles 
of the same quantities in the CO simulation at f = 250 Myr and 
t = 500 Myr Squares have been added to the profiles of the CO 
simulation at f = 250 Myr in order to aid comparison with the 
C02 case. 



of the cluster due to cooling. Profiles are obtained for the C02 
(solid lines) and CO (dashed lines) simulations by averaging on 
density and assuming a bin interval of 3 kpc. The temporal evo- 
lution for the C02 simulation is represented by the line thick- 
ness, which decreases with time: profiles are computed at inter- 
vals of 50 Myr from the initial conditions up to f = 250 Myr 
The dashed lines represent the profiles of the same quantities in 
the CO simulation to t - 250 Myr and t - 500 Myr; at t - 250 
Myr we also added squares to the dashed line to facilitate visual 
comparison. The figure shows that the profiles at / = 250 Myr 
are practically coincident to the limit of resolution; this implies 
that the finite size of the box does not induce spurious effects in 
the C02 simulation. 

As expected, these runs show the establishment of an ho- 
mogeneous cooling flo w which evolves in a runaway manner 
dWhite & Sarazinl 11987 ). Because of the long cooling time, it 
is not surprising that the temperature of the core drops only to 
r = 3 X 10^ K after the first 250 Myr Using a 2D simulation 
with the same resolution as the C02 case, we checked that the 
central temperature falls to the minimum value between t - 390 
and 400 Myr. We will use the profiles obtained in the C02 case 
as references for the same quantities computed in the presence 
of bubbles. 

3.2. The rise of the bubbles 

The rise of spherical bubbles of plasma in a cluster atmosphere 
has been described in many papers (e.g. Churazov et al., 2000). 
We observe that the evolution of bubbles in the ICM resem- 
bles quite closely the dynamics of powerful explosions in the 



earth's atmosphere, or more generally of "thermals", which are 
spherical bubbles of hot air in pressure equilibrium with the sur- 
rounding colder air In the latter case, the main difference resides 
in the fact that thermals are only slightly hotter and less dense 
than the surrounding air; this produces a rather different behav- 
ior, and thermals are able to rise in the atmosphere much higher 
than their initial diameter, expanding linearly with time ( Turneii 
1973). Studies of thermals can hence adopt the Boussinesq ap- 
proximation in computations and perform experiments using liq- 
uids. Nevertheless, some of the behavior of thermals, such as 
their internal motions, presents remarkable similarities with the 
behavior of plasma bubbles. The behavior of spherical bubbles 
in water also resembles, for a p roper choice of parameters, the 
behavior of plasma bubbles (e.g. [Walters & Davidsonl[l963h . 

The evolution of the 3D simulations with bubbles is shown 
in Figures 4 to 9 through maps constructed by cutting the sim- 
ulations along one of the planes of symmetry which contain the 
centers of the cluster and the bubbles. We concentrate on the 
simulation SB 1011 as representative of the sample, and will re- 
fer to the other 3D simulations only when their behavior is quite 
different with respect to it. 

We paid some attention to the first stages of the evolution, 
namely the starting of the bubble from rest and the generation of 
sound waves. This phase is usually n eglected in the literatur e as 
transient and unrealistic (see however [Churazov et al.L 1200 2*) de- 
spite the fact that sound waves are often clearly visible in figures. 
In Figure 4 we show the evolution of the pressure and velocity 
fields during the first 10 Myr The plots show how the bubble 
plasma oscillates vertically and generates sound waves, while 
motions are driven into the adjacent ICM at the same time. At 
the time of the last snapshot, a vortex circulation directed up- 
ward along the bubble axis a nd downward on the bubble surface 
dffiiil 1 1 894t iBatchelod. [T967h is established in the plasma and in 
the surrounding ICM. This kin d of circulation enables the bubble 
to rise steadily ( lTurneitri973h . 

During the rise that follows the initial transient phenomena, 
each bubble is affected by instabilities on its surface due to lack 
of surface tension. In particular it will suffer Rayleigh-Taylor 
(RT) instabilities on the top and Kelvin-Helmoltz (KH) instabil- 
ities on its flanks. In the end, however, it is the vortex circulation 
itself which mostly affects the bubble. Indeed, in an ideal case 
a stagnation layer should be established in the ICM below the 
bubble, but in reality, as the bubble rises, the fall of pressure at 
its bottom draws the ICM upward, where it is entrained in the 
vortex circulation, thus acting as an effective RT instability. 

The evolution of the rising bubbles is depicted in Figure 5 
with intervals of 10 Myr up to f = 120 Myr and in Figures 6 to 
8 with intervals of 50 Myr up to f = 250 Myr. As long as the 
bubble rises, a trunk of denser material is entrained in its wake, 
which penetrates it from below and gives it a mushroom shape. 
As expected in the development of RT instabilies, the uplifted 
material rises faster and eventually bisects the bubble, which 
then assumes the shape of a rotating torus (or vortex ring). It 
is worth noting that the ICM entrained by the bubble cools by 
expansion because of the lower pressure, and at some times it 
is the coldest material in the cluster. As the evolution goes on, 
the lower part of the trunk falls back to the cluster center, while 
the upper part remains at large radii as long as it is driven in the 
vortex circulation. 

Defining the borders of the bubbles becomes a harder task 
as the simulation proceeds. Figure 8 shows the fraction of the 
bubble plasma in the fluid elements. Inspection of these values 
shows clearly that all of the hot plasma has mixed completely 
with the ICM after 100 Myr In particular, the plasma picked up 




Fig. 4. The evolution of the pressure field from f = Ouptof = 10 Myr for the simulation SB 101 1 in the region around the bubble, 
with the velocity field superimposed. The first picture shows the border of the bubble as well as a velocity vector corresponding to 
V = 1000 km/s. The last snapshot shows the region occupied by the fluid elements containing at least 90% hot plasma at f = 10 Myr, 
with the initial position of the bubble superimposed on it. The bubble rises ~ 2 kpc in the first 10 Myr. 



along the bubble axis by the rising ICM cools rapidly to low tem- 
perature due to the high density of the entrained ICM and should 
no longer be treated as relativistic. In the vortex ring, the plasma 
is also mixed with the ICM, but because of the lower density of 
the polluting ICM, the mixture still preserves buoyancy and a 
high temperature. 

As stated above, despite its relevance for the evolution of the 
bubbles and the ICM dynamics, mixing is not properly treated 
in these simulations because the diffusion of species is not ex- 
plicitly accounted for. We address this point by checking that 
the amount of mixing remains stable with respect to increases in 



resolution. We ran a 2D simulation with the same characteristics 
as the SB1020 case and compared it with an identical 2D sim- 
ulation with two additional levels of refinement. We found that 
both simulations show the same amount of mixing up to f = 250 
Myr, despite the fact that the plasma distribution shows narrower 
features in the second case, as expected. We consider that the 
mixing should eventually increase in the 2D simulations due to 
their enhanced sensitivity to instabilities; hence the fact that the 
amount of mixing remains the same is a proof of its robustness. 

The vortex rings evolve similarly in the different simulations. 
In all of the models the rings are affected by instabilities at reg- 




Fig. 5. Map of the temperature in the SB 1011 simulation showing the evolution of the plasma bubble in the first 120 Myr. The 
pictures show clearly the effect of the entrained cold gas, which penetrates and finally bisects the ICM bubble, leaving a vortex ring. 



ular distances along their length that act to break the rings into 
smaller sections. It is worth noting that the positions of the break 
points along the rings are symmetric with respect to the princi- 
pal vertical planes of the system, and hence the instabilities do 
not reflect the axial symmetry of the cluster but the geometry of 
the computational grid. In the SB 1011 simulation, the rings are 
already broken in the last outputs exactly along the plane of the 
maps, so that it seems in Figure 6 that only some remnants are 
slightly hotter than the surrounding medium after t - 200 Myr 
However, this is not the case: the hot material is concentrated 
outside of the plane of the figure. In the SB1020 and SB1030 
simulations, the vortex rings are not yet broken by the end of 
the run, but instabilities are effectively breaking them outside 



the plane of the image. In all of the simulations the rings pre- 
serve high temperatures despite the fact that the fraction of bub- 
ble plasma inside them is less than 10%. 

We quantify the evolution of the bubbles in Figure 10, which 
shows the height of the centroid of the plasma material as a func- 
tion of time. The dashed lines describe the value of the height of 
the centroid of all of the plasma, while the solid lines depict the 
same quantity for the plasma that is hotter than T = 8 x 10'' K 
(the maximum temperature of the ICM inside r = 80 kpc). The 
thickness of the lines distinguishes among the different simula- 
tions. The three simulations show remarkably similar behavior: 
the hotter material, which effectively defines the bubbles, rises 
by ~ 20 kpc in the first 100 Myr, and then decreases its speed 
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Fig. 6. The evolution of the ICM temperature in the SB 1011 
simulation from the beginning of the simulation up to f = 250 
Myr The bubble plasma seems almost completely diffused into 
the ICM in last snapshot, but it is still preserved in remnants of 
the vortex ring that do not intersect the plane of the image. The 
pictures show the entrainment of cold gas along the vertical axis, 
the establishment of a cold cap at large radii, the fall back of the 
cold gas and the development of the cooling flow. 



until the end of the simulation. The average speed is then ~ 200 
km/s for the first 100 Myr, which can be compared to the speed 
of sound in the ICM, ranging from ~ 1000 km/s at the center to 
~ 1500 km/s at r = 100 kpc. The differences between the values 
for the overall plasma and the hottest plasma at the beginning of 
the simulation are explained by the fact that, due to the finite grid 
size, a portion of bubble material is mixed with the ICM akeady 
in the initial conditions. Initially this material is less aff'ected by 
buoyancy, and it is driven upward later in the bubble wake. The 
rising of the trunk of cold gas, the sweeping of the plasma on 
the bubble axis and the establishment of a cold cap above the 
ring, are also well represented by the rise of the overall centroid 
higher than the hot plasma. The succeeding decrease depicts the 
falling back of the cold material along the trunk. 

The centroid of the hot material still seems to rise steadily 
in the last outputs of our simulations, but inspection of the frac- 
tion of plasma that is still hot, depicted in Figure 11, shows that 
the plasma in the SB 1011 simulation is rapidly cooling. In the 
SB 1020 and SB 1030 simulations, it instead appears that the frac- 
tion of hot plasma is stable around 35% in the final output. 



Fig. 7. The evolution of the density in the SB 1011 simulation 
up to f = 250 Myr. The snapshots show the denser material that 
is uplifted by the vortex circulation and later sustained at large 
radii by the buoyancy of the vortex ring. 

To explore the later stages of the evolution, we turn again to 
2D simulations with the same parameters and resolution as the 
3D cases. In general, these simulations show that material con- 
tinues to fall along the trunks down to the cluster center, while 
the ICM already present there is displaced orthogonally and set- 
tles into a torus. The trend then reverses and the material in the 
torus falls back onto the core, driving back upward the mate- 
rial that has fallen along the trunks. Overall, this resembles a 
sort of oscillatory motion, but it becomes more turbulent as the 
simulation proceeds. We halt the simulations when some plasma 
element cools close to the minimum temperature and the com- 
putational time step becomes too short; we can suppose that a 
runaway cooling similar to the homogeneus cooling flow is es- 
tablished at that time. Vortex rings are preserved until the end of 
the run, deforming and slowly rising for few kpc more. 

4. Effects of bubbles on the cooling 

It has been observed that AGNs might heat the ICM and pos- 
sibly quench cooling flows in any of several different ways. 
In particular, it has bee n predicted (e.g. Heinz et al., 1998) and 
shown by simulations (iZanni et al.L 12005: .Briiggen et aU ,20051) 
that most of the jets' energy is transferred to the ICM through 
the bow shock and the sound waves pushed by the expanding 
cocoon. These processes deposit energy into the ICM as the jets 
cross it, affecting the ICM at a distance from the radio source 
itself and possibly outside the cool core, unless some viscosity 
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Fig. 8. The fraction of bubble plasma in the fluid elements in 
the SB 1011 simulation up to f = 250 Myr. The comparison with 
Figure 6 and 7 shows that as the simulation goes on, most of 
the bubble plasma is indeed mixed with the ICM to low temper- 
ature and high density, and under these conditions it should no 
longer be treated as relativistic. The reader should be careful to 
the logarithmic colour scale. 



Fig. 9. 

tion. 



The evolution of the temperature in the SB 1020 simula- 



or thermal conductivity is present ("Ruszk owski et al.l l2004allbl : 
[Briiggen et al., 2005; Fabian et al., 2005). On the other hand, en- 
ergy is still stored in the cavities as thermal energy and in the 
ICM in potential form during the later phase of evolution by 
buoyancy. Both these forms of energy are released during the 
rise of the cavities: the first through pdV work performed by 
the cavities, which expand in order to reach or mainta i n pres - 
sure equilibrium with the environment (e.g. iBegelmanl 1200 1[) . 
and the second when the surrounding medium falls in around 
the cavity to fil l the space the bubb l e prev iously occupied (e.g. 
IChurazov et"an. l2002t iBirzan et all [2004). Work and thermal 
energy exchange also occur when portions of the ICM at differ- 
ent temperatures and pressures come into contact and mix, after 
being displaced by motion s induced by the rise of the cavities 
(Quilis e t al, 2001; Churaz ov etaTl 120021: iDaUa Vecchia et all 
|2004|) . This mechanism becomes more eff'ective in heating the 
cool core as more material is radially displaced. 

Despite the presence of bubbles, our simulations show that 
a cooling flow is established in the central regions of the cluster 
(Figures 6 and 9). This is even more true as the initial location 
of the bubbles is moved away from the center This is hardly 
surprising, given that bubbles are the outcome of a single burst 



of energy injection while cooling is continously operating in the 
ICM. Nonetheless, bubbles affect the ICM and can slow down 
its cooling effectively; in the following we investigate how and 
by what amount. 

Because of the dependence of gas emissivity upon density 
(approximately e oc p^T""^), cooling is enhanced by deposition 
and compression of material. Thus, we observe that bubbles af- 
fect the cooling firstly by their mere presence, because of the 
removal from their volumes of some fraction of the cold ICM 
akeady stored in the core. Even if this can be regarded as an ar- 
tifact of our setup, we expect that a similar effect is produced on 
the ICM by displacement and shock heating during bubble infla- 
tion. Later, some amount of thermal energy is directly deposited 
in the ICM through mixing with the bubble material, while stir- 
ring of the ICM due to the bubble rise is expected to reduce the 
deposition of cooling material to the core. To investigate the ef- 
fects of the gas motions on the cooling, we performed a new run 
of the SB 1020 case using a tracer fluid to mark the fluid elements 
initially stored inside r = 10 kpc; we show in Figure 12 the evo- 
lution of their distribution. We observe that in the first stages 
of evolution, the central material is stretched and uplifted in the 
bubbles' wake and expands rather than being compressed. Later 
on, material indeed falls onto the core, but because of the geom- 
etry, it accumulates directly at the center rather than on the cold 
ICM, which is displaced instead into a torus. A similar analysis 
was performed also in the SB 101 1 case; its results are similar to 
SB 1020 up to r = 190 Myr, when the run halts. 
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Fig. 10. The height of the centroid of the plasma material during 
each run for the complete set of bubble simulations. The dashed 
Unes describe the height of the centroid of all the bubble plasma, 
while the solid lines depict the same quantity for the plasma that 
is hotter than T - 8 X 10^ K, which is the maximum temperature 
inside r = 80 kpc. 

We measure the effects of the bubbles on the cooling by 
comparing the temperature profiles with the profiles of the C02 
model. This is done in Figure 13, which reports the difference 
T - Tco2 from the initial conditions up to f = 250 Myr with inter- 
vals of 50 Myr. It must be remembered that, as depicted in Figure 
3, in the absence of bubbles the gas cools about AT ~ 2 x 10^ K 
during the first 250 Myr of evolution, while the cooling is much 
faster at later times. Hence, a slight gain in the core temperature 
can be indicative of a significantly longer evolution of the cool- 
ing process. The bubbles in simulation SB 1011 have the great- 
est effect on the cooling, as expected, due to their proximity to 
the core of the cluster. In this case the temperature is increased 
with respect to the C02 simulation in all of the region inside 
r ~ 50 kpc for most of the run. Att = 250 Myr, the tempera- 
ture difference reaches T - Tcoi ~ 10^ K. The most interesting 
result comes from the SB 1020 simulation. In this case the cool- 
ing in the center proceeds apparently as in the C02 case up to 
t ~ 200 Myr. After this time, however, the central temperature 
drops more slowly, and in the final snapshot T - Tcoi ~ 6 x 10^ 
K. This is due to the replacement of the cold material in the core 
by the dense gas falling to the center along the trunks. In the 
SB 1030 simulation the core of the cluster instead appears com- 
pletely unaffected by the presence of bubbles up to the end of the 
run at f = 250 Myr. However, because dense gas is clearly falling 
along the trunks at that time, a later heating as in SB 1020 can- 
not be excluded. Because of the fraction of plasma entrained in 



Fig. 11. The fraction of bubble plasma which maintains a tem- 
perature r>8xlO^Kasa function of time. Even if the tem- 
perature of this plasma is lowered to the non-relativistic regime, 
it is stiU hot and underdense enough to be driven by buoyancy. 

the core (Fig. 8), we wondered if the higher temperature of the 
infalling material was due to mixing. Therefore we computed 
the profiles of the fraction of internal energy due to the bubble 
plasma as shown in Figure 14. At ? = 250 Myr, we found that the 
internal energy directly deposited in the ICM amounts at most at 
2% everywhere except in the shells hosting the vortex rings. 

The use of 2D simulations enables us to extend this analy- 
sis up to the establishment of runaway cooling when some gas 
element reaches the mimimum temperature. In the case of the 
homogeneus cooling flow, we saw that this happens soon after 
/ - 390 Myr. In SB 1011, after a turbulent oscillation, runaway 
cooling is established in the core at some time after t = 500 Myr. 
SB1020 presents an analogous behavior, but because of the later 
infall, the cooling flow is more rapidly re-established and the 
runaway phase starts at f = 440 Myr. In the SB 1030 case, the in- 
falling material does not reach the center but rather compresses 
the material already present in the core; the effect is to enhance 
the cooling so that the runaway phase begins after ? = 380 Myr. 

5. Conclusions 

We show the results of a set of 3D simulations of buoyant bub- 
bles in the atmosphere of a cool core cluster. We allow cooling, 
consider the bubble plasma as relativistic, and achieve high ef- 
fective resolution. We consider the presented simulations to be a 
first step in studying the interaction of plasma bubbles with the 
ICM. We will use these runs as references for future work with 
different configurations of the system and additional physics. 
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Fig. 13. The difference between the temperature profile of the 
simulations hosting bubbles and the C02 case. 



Fig. 12. The evolution of the material initially inside r = 10 
kpc in run SB 1020. It is first stretched by the rise of the bubbles 
and later displaced into a torus around the center by the material 
falling along the trunks. 



The d ynamics of our bub bles res embles the res u hs de - 
scribed by ISaxton et al.1 (l200ll) and by [Reynolds etaTI (l2005l) . 
ISaxton et al.l (1200 ll ) produce a set of 2D simulations similar to 
SB 1030 for difi'erent values of the ratio of the density of plasma 
versus ICM. They remark upon the robustness of the vortex rings 
and the fall of entrained material to the core along the trunks; as 



in SB 1030, this material simply depos its onto the core rather 
than disrupts it. ' Reynolds et al.l ( 1200 5^ study the evolution of 
bubbles for different values of the ICM viscosity using 3D simu- 
lations. In a reference case without viscosity the bubble evolves 
similarly to SB 101 1, deforming in vortex rings with a thin trunk 
of material entrained in its wake. Each of these studies focuses 
on the evolution of bubbles, rather than on their impact on the 
environment, and therefore neglects cooling. 

Our simulations do not fully address the problem of the 
quenching of the cooling flow by AGN, because the phase of ac- 
tivity is neglected and many simplifications are assumed in the 
setup. We address, however, the direct effect of our bubbles on 
the cooling and find that it is rather modest. Even in the SB 101 1 
case, runaway cooling is postponed only by about 25% respect to 
the unperturbed case. Setting the bubbles further away only re- 
duces their effect, or even accelerates the cooling as in SB 1030. 
Larger bubbles will probably have a big ger impact, but the di- 
mensions chosen in our cases are typical dBirzan et ani2004l) . 



Numerically, the results are quite robust with respect to in- 
crease in resolution and the adoption of different systems of 
coordinates. Thus we conclude that our resolution is adequate. 
Details are different, however, and 2D simulations are much 
more sensitive to the propagation of rounding errors. For this 
reason we prefer to trust the 3D cases and use 2D simula- 
tions only to extend their results. To a lesser extent, 3D sim- 
ulations suffer numerical instabilities as well, as revealed by 
the shredding of vortex rings along the lines of the computa- 
tional grid. One technique for mitigating grid effects involves 
just adding some random noise in the initial conditions, as in 
i Jones & De Younj (l2005h . The symmetry of our results, how- 
ever, testifies to the reliability of the code. 

Future work will consider some relevant physics that is not 
included in the present simulations. In particular, we will con- 
sider the magnetic field and the ICM viscosity, which have been 
shown to be able, for proper values, to stabilize the bubble sur- 
face. Using more realistic initial conditions, we will take into 
better account the activity phase of AGN and its effects on the 
medium. The introduction of thermal conduction can be relevant 
in heating the entrained cold gas as it rises to regions of higher 
temperature, besides directly heating the cluster core. Mixing 
of hot and cold gas should also be properly treated, as well 
as the transition of the plasma to the non-relativistic regime. 
We aim to study also the interaction of pre-existing bubbles 
with shock waves induced by AGN activity or merger events 
(EnBIin&Bruggen, 2002; Heinz _& Churazov, 2005). Finally, 
a direct comparis on with X-ray observations is possible using 
tools like X-MAS dGardini etalil2004 . 
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